COMPUTATION OF EXTERIOR MODULI OF QUADRILATERALS 



HARRI HAKULA, ANTTI RASILA, AND MATTI VUORINEN 

Abstract. We study the problem of computing the exterior modulus of a bounded 
quadrilateral. We reduce this problem to the numerical solution of the Dirichlet- 
Neumann problem for the Laplace equation. An algorithm and its implementation 
with ft,p-FEM is documented. Several experimental results, with error estimates, 
are reported. 



1. Introduction 

A bounded Jordan curve in the complex plane divides the extended complex 
plane Coo = C U {00} into two domains Di and D2, whose common boundary it 
is. One of these domains, say Di , is bounded and the other one is unbounded. The 
domain Di together with four distinct points zi,Z2,Z3, Z4 in dDi , which occur in this 
order when traversing the boundary in the positive direction, is called a quadrilateral 
and denoted by {Di, zi,Z2,Z3, Z4) . 

By Riemann's mapping theorem, the domain Di can be mapped conformally 
onto a rectangle f : Di ^ (0, 1) x (0, h) such that the four distinguished points are 
mapped onto the vertices of the rectangle f{zi) = , f{z2) = 1, fiz^) = 1 + ih, 
f{z4) = ih . The unique number h is called the (conformal) modulus of the quadrilat- 
eral [Di] zi, Z2, z^, Z4) . Apart from its theoretical significance in geometric function 
theory, the conformal modulus is closely related to certain physical quantities which 
also occur in engineering applications. In particular, the conformal modulus plays 
an important role in determining resistance values of integrated circuit networks (see 
e.g. [221 ES]). Similarly, one can map D2 , the complementary domain, conformally 
g: D2 ^ (0; 1) ^ i^jk) such that the four boundary points are mapped onto the 
vertices of the rectangle g{zi) = 0, g{z2) = 1, g^z^) = 1 + ik,g{z4) = ik , reversing 
the orientation. Again the number k is unique and it is called the exterior modulus of 
{Di, zi, Z2, zs, Z4) . In practice, the computation of both the modulus and the exterior 
modulus is carried out by using numerical methods such as numerical conformal map- 
ping. Mapping problems involving unbounded domains likewise are related to some 
well known engineering applications such as determining two dimensional potential 
flow around a cylinder, or an airfoil. 

In the case of domains with polygonal boundary, numerical methods based on 
the Schwarz-Christoffel formula have been extensively studied, see [TO] • The literature 
and software dealing with numerical conformal mapping problems is very wide, see 
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e.g. [in] and [22]. In our earlier paper [H] we applied an alternative approach which 
reduces the problem to the Dirichlet-Neumann problem for the Laplace equation. 
Thus any software capable of solving this problem may be used. We use the hp-FEM 
method for computing the modulus of a bounded quadrilateral and here we will 
apply the same method for the exterior modulus and another method, AFEM, for 
the sake of comparison, as in [H]. Our approach also applies to the case of domains 
bounded by circular arc boundaries as we will see below. It should be noted that 
while our method does not require finding the canonical conformal mapping, it is 
possible to construct the mapping from the potential function. An algorithm, with 
several numerical examples, is presented in [13\ . 

In particular, an important example of a quadrilateral [Di] zi, Z2, z^, Z4) is the 
case when Di is a polygon with ^1,^2, ^3, ^4 as the vertices and its modulus was 
computed in [Kj and this formula was also applied in [13] . Here we reduce its exterior 
modulus to the (interior) modulus by carrying out a suitable inversion which keep 
three vertices invariant and maps the exterior to the interior of a bounded plane 
region whose boundary consists of four circular arcs. 

We apply here three methods to study our basic problem: 

(1) The hp-FEM method introduced in fT^ and its implementation by H. Hakula. 

(2) The AFEM method of K. Samuelsson, see e.g. [7] and |14j . 

(3) The Schwarz-Christoffel Toolbox of T. DriscoU and N. Trefethen P [TOl EHl ET] . 

The methods (1) and (2) are based on a reduction of the exterior modulus problem 
to the solution of the Dirichlet-Neumann problem for the Laplace equation in the 
same way as in [T] and [E] whereas (3) makes use of numerical conformal mapping 
methods. Note that [T] also provides a connection between the extremal length of a 
family of curves and its reciprocal, the modulus of a curve family, both of which are 
widely used in the geometric function theory. 

We describe the high-order p-, and /ip-finite element methods and report the 
results of numerical computation of the exterior moduli of a number of quadrilat- 
erals. In the p-method, the unknowns are coefficients of some polynomials that are 
associated with topological entities of elements, nodes, sides, and the interior. Thus, 
in addition to increasing accuracy through refining the mesh, we have an additional 
refinement parameter, the polynomial degree p. For an overview of the hp-method, 
see e.g. Babuska and Suri p]. A more detailed exposition of the methods is given in 

[23 Eg. 

Our study is structured according to a few particular cases. We start out with 
the case when the quadrilateral is the complement of a rectangle and the vertices 
are the distinguished points of the quadrilateral. In this case we have the formula of 
P. Duren and J. Pfaltzgraff [llj to which we compare the accuracy of each of the above 
methods (l)-(3). Then we consider the problem of minimizing the exterior modulus 
of a trapezoid with a fixed height h and fixed lengths for the pair of parallel opposite 
sides and present a conjecture supported by our experiments. We also remark that the 
case of symmetric hexagons can be dealt with the Schwarz-Christoffel transformation 
and relate its exterior modulus to a symmetry property of the modulus of a curve 
family. Finally, we study the general case and present comparisons of (l)-(3) for 
this case as well. SC Toolbox does not have a built in function for computing the 
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exterior modulus. However, we use the function extermap, and an auxiliary Mobius 
transformation, to map the exterior of a quadrilateral (Z); a, b, c, d) conformally onto 
the upper half-plane so that the boundary points a, b, c and d are mapped to the points 
oo, —1,0 and t > 0, respectively. Then the exterior modulus of the quadrilateral is 
r(t)/2, where r is the Teichmiiller modulus function (see ^). We use the MATLAB 
code from [2] to compute values of r(t), t > 0. 

At the end of the paper we present conclusions concerning the performance of 
these methods and our discoveries. 

2. Preliminaries 

In this section we give reference results which can be used in obtaining error 
estimates. We also present some geometric identities which are required in our com- 
putations. 

2.1. The hypergeometric function and complete elliptic integrals. Given 
complex numbers a,b, and c with c 7^ 0,-1,-2,..., the Gaussian hypergeometric 
function is the analytic continuation to the slit plane C \ [1, 00) of the series 

(2.2) F(a,6;c;z)=2Fi(a,6;c;z) = f;^^^;^^^, \z\ < 1 . 

n=0 ^ ' ' 

Here (a, 0) = 1 for a 7^ 0, and (a, n) is the shifted factorial function or the Appell 
symbol 

(a, n) = a{a + l)(a + 2) ■ ■ ■ (a + n — 1) 

for G N \ {0}, where N = {0, 1, 2, . . .} and the elliptic integrals X(r), %'{r) of the 
first kind are 

X(r) = 1^(1/2, 1/2; 1; r^), X'{r) = X(r'), and r' = Vl - r^, 
and the elliptic integrals £(r), £'(r) of the second kind are 
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£(r) = 2 ^(1/2, -1/2; 1; r^), e'(r) = £(r'), and r' = Vl - r 
Some basic properties of these functions can be found in [21 [21] . 

2.3. The modulus of a curve family. For a curve family F in the plane, we 
use the notation M(F) for its modulus [20] • For instance, if F is the family of all 
curves joining the opposite 6-sides within the rectangle [0, a] x [0,6], a, 6 > 0, then 
M(F) = b/a. If we consider the rectangle as a quadrilateral Q with distinguished 
points a + ib, ib,0,a we also have M (Q; a + ib, ib,0,a) = b/a , see [H I2Q] . Given three 
sets D,E,F we use the notation A[E, F; D) for the family of all curves joining E 
with F in D . 

2.4. The Duren-Pfaltzgraff formula [TU Theorem 5]. For k G (0, 1) write 

_ 2(£(fc)-(l-A:)X(fc)) 
^ £'(A;) - kX'ik) • 

Then ip: (0, 1) — t- (0, 00) defines an increasing homeomorphism with limiting values 
0,00 at 0,1, respectively. In particular, ip~^ : (0,oo) — )■ (0,1) is well-defined. Let R 
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be a rectangle with sides of lengths a and b, respectively, and let T be the family of 
curves lying outside R and joining the opposite sides of length b . Then 



(2.5) 



M(r) 



X'{k) 



2%{k) ' 

See also W.G. Bickley P (1.17) p. 86]. 



where k = ip '^{a/b) . 



2.6. Circle through three points. First we find the equation of the circle {x — 
^oY + {y ~ yoY = ^"^ through the points Zj = Xj + ii/j, j = 1, 2, 3. It is well known 
that 
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a,b,0 onto the unit 



Now we may define a similarity mapping which maps zi,Z2, z^ 
circle by 

w{z) = {z- ZQ)/r, 
where Zq = Xq + iyo and r is as above. 

2.7. Mapping unbounded onto bounded domains. The transformation z i— )■ 
z/lzl"^ maps the complement of the closed unit disk onto the unit disk. This trans- 
formation is an anticonformal mapping and it maps the complement of a polygonal 
quadrilateral with vertices a,b,c,d with |6| = |c| = = 1 onto a bounded domain, 
bounded by four circular arcs. Note that the points b, c, d remain invariant under 
this transformation. See Figure [T} 

2.8. The Dirichlet-Neumann problem. The following problem is known as the 
Dirichlet- Neumann problem. Let Dhe a. region in the complex plane whose boundary 
dD consists of a finite number of regular Jordan curves, so that at every point, except 
possibly at finitely many points, of the boundary a normal is defined. Let dD = AUB 
where A,B both are unions of Jordan arcs. Let ipA^i^B be a real-valued continuous 
functions defined on A,B, respectively. Find a function u satisfying the following 
conditions: 

(1) u is continuous and differentiable in D. 

(2) u{t) = ipAit), for all t e A. 

(3) If d/dn denotes differentiation in the direction of the exterior normal, then 

d 

—uit) = iPbU), for all t e B. 
on 
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Figure 1. Polygonal quadrilateral before (left) and after (right) the 
inversion transformation z ^ zj\z\^ . Note that the points 6, c, d on the 
unit circle remain invariant. 

2.9. Modulus of a quadrilateral and Dirichlet integrals. One can express the 
modulus of a quadrilateral (D\ zi, Z2, z^, Z4) in terms of the solution of the Dirichlet- 
Neumann problem as follows. Let 7^, j = 1, 2, 3, 4 be the arcs of dD between (^4, zi) , 
{zi,Z2) , (-2^2, -23) , (-23,^4), respectively. If u is the (unique) harmonic solution of the 
Dirichlet-Neumann problem with boundary values of u equal to on 72, equal to 1 
on 74 and with 9M/(9n = on 71 U 73 , then by [U p. 65/Thm 4.5]: 



2.11. The reciprocal identity. Given a quadrilateral Q = {D; zi,Z2, ^3, Z4) we call 

sometimes Q = {D; Z2, Z3, z^, Zi) the conjugate quadrilateral. It a simple basic fact 
that 



might serve as a useful error characteristic. We will continue to use this also in our 
work. 



In this section we describe the high-order p-, and hp-Unite element methods. The 
paper of Babuska and Suri [B] gives an accessible overview of the method. For a more 
detailed exposition we refer to Schwab [23], and for those familiar with engineering 
approach the book by Szabo and Babuska [2S] is recommended. The presentation 
here follows closely earlier work by the authors [T4] . 

For the applications considered in this paper, any finite element computation 
requires at least the choice of the following. 



(2.10) 





3. The hp-FEM method 
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(1) Initial discretization of the domain. In 2D each discretization or mesh divides 
the domain into elements, plane regions with piecewise smooth boundaries. 
These are usually either triangles or quadrilaterals. 

(2) Refinement strategy. The choice of the refinement strategy is connected to 
choosing the finite element method (FEM): mesh refinement (/i-method), ele- 
mentwise polynomial order (p-method), or both (hp-method). The unknowns 
or degrees of freedom are the coefficients of the chosen shape functions. In the 
/z.-version, the shape functions are such that the coefficients are also values 
of the solution at specified locations of the discretization of the computa- 
tional domain, that is, the nodes of the mesh. In the p-method, the shape 
functions are polynomials that are associated with topological entities of the 
elements, nodes, sides, and interior. Thus, in addition to increasing accuracy 
through refining the mesh, we have an additional refinement parameter, the 
polynomial degree p. 

Refinement strategies can be based on a priori and a posteriori error estimates. 
In all cases considered here, sufficient a priori knowledge of the singularities exist for 
the construction of the properly graded meshes. All modifications of the refinement 
strategy always imply recomputing of the entire solution. 

What makes the exterior modulus problem especially interesting from the finite 
element point of view is the natural appearance of computational domains with cusps. 
Our experience shows that the standard geometric grading is adequate also in this 
class of problems. In the cases when reference solutions are not available, we rely on 
the reciprocal identity as our fundamental quality measure. 

Let us next define a p-type quadrilateral element. The construction of triangles 
is similar and can be found from the references given above. 

3.1. Shape functions. Many different selections of shape functions are possible. We 
follow Szabo and Babuska [25] and present the so-called hierarchic shape functions. 
Legendre polynomials of degree n can be defined using a recursion formula 



For our purposes the central polynomials are the integrated Legendre polyno- 
mials for X G [—1, 1], 




(3.3) 




n = 2,3, . . . , 



which can be rewritten as linear combinations of Legendre polynomials 



(3.4) 0„(^) = ^^,^(P„(^)_P„_2(^)), n = 2,3,.... 




The normalizing coefficients are chosen so that 



(3.5) 
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We can now define the shape functions for a quadrilateral reference element over 
the domain [—1, 1] x [—1, 1]. The shape functions are divided into three categories: 
nodal shape functions, side modes, and internal modes. 

3.2. Nodal shape functions. There are four nodal shape functions: 













^(1 + 0(1 + 7?), 




^(1-0(1 + ^)- 



Taken alone, these shapes define the standard four-node quadrilateral finite element. 

3.3. Side shape functions. There are 4(p — 1) side modes associated with the sides 
of a quadrilateral (p > 2) : 







1 = 2,.. 


■ ■ ,P, 






i = 2,.. 


■■,P, 






1 = 2,. 


■■,P, 






1 = 2,.. 
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3.4. Internal shape functions. For the internal modes we have two options. The 
so-called trunk space has (p — 2)(p — 3)/2 shapes 

(3.6) Nl^i^,r]) = UOMv), i,j>'^, i + J = 4, 5, . . . ,p, 
whereas the full space has (p — l)(p — 1) shapes 

(3.7) Nl^{^,n) = MOMv): ^ = 2,...,p, i = 2,...,p. 

In this paper we always use the full space. The internal shape functions are often 
referred to as bubble-functions. 

3.5. Parity problem. The Lcgcndrc polynomials have the property Pn{—x) = 
(— l)"Pji(a;). In 2D all internal edges of the mesh are shared by two different ele- 
ments. We must ensure that each edge has the same global parametrization in both 
elements. This additional book-keeping is not necessary in the standard /i-FEM. 
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Figure 2. Curved boundary mapping. 



3.6. Domains with curved boundaries. It is important to represent curved bound- 
ary segments accurately. Tlie linear blending function method of Gordon and HaU 
[12] is our choice for this purpose. 

In the general case all sides of an element can be curved as in Figure [2j We 
assume that every side is parametrized as follows: 

(3.8) X = Xi{t), y = Uiit), — 1 < t < 1, i = I, . . . , number of sides 

Using capital letters as coordinates of the corner points, {Xi,Yi), we can write the 
mapping for the global x-coordinates of a quadrilateral as 

^ = ^(1 - v)MO + ^(1 + O^^iv) + ^(1 + v>3iO + ^(1 - OMv) 

(3.9) - ^(1 - 0(1 - v)Xi - ^(1 + 0(1 + V)X2 - ^(1 + 0(1 + v)Xs 
-^(l-0(l + r])X4, 

and symmetrically for the y-coordinate. Note, that if the side parametrizations 
represent straight edges, the mapping simplifies to the standard bilinear mapping of 
quadrilaterals. 

In the following we always use exact representation of the geometry which im- 
plies that in the ensuing mesh grading process no approximation of geometry is 
necessary. 

3.7. Proper grading of the meshes. For a certain class of problems it can be 
shown that if the mesh and the elemental degrees have been set optimally, we can 
obtain exponential convergence. A geometric mesh is such that each successive layer 
of elements changes in size with some geometric scaling factor a, toward some point 
of interest. In this case, the points of interest are always corner points. 

The following theorem is due to Babuska and Guo [Ij. Note that construction 
of appropriate spaces is technical. For rigorous treatment of the theory involved see 
Schwab [21], Babuska and Guo [5J and references therein. 

Theorem 1. Let ^7 C be a polygon, v the FEM-solution, and let the weak solution 
uq be in a suitable countably normed space where the derivatives of arbitrarily high 
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Figure 3. A sample geometry and the corresponding initial mesh. 
Note that the three-element -rule is satisfied at every corner. 

order are controlled. Then 

inf \\uo -vWh^q) < C exp(-6v^), 

V 

where C and b are independent of A^, the number of degrees of freedom. Here v 
is computed on a proper geometric mesh, where the orders of individual elements 
depend on their originating layer, such that the highest layers have the smallest 
orders. 

The result also holds for constant polynomial degree distribution. 

Let us denote the number of the highest layer with u, the nesting level. Using 
this notation we refer to geometric meshes as {a, z/)-meshes. 

3.8. Generating geometric meshes. We use the following two-phase algorithm 
where triangles can be replaced by quadrilaterals or a mixture of both: 

(1) Generate an initial mesh (triangulation) where the corners are isolated with 
a fixed number of triangles depending on the interior angle, 6 so that the 
refinements can be carried out independently: 

(a) 9 < 7i/2: one triangle, 

(b) 7i/2 < 9 < Ti: two triangles, and 

(c) TT < 9: three triangles. 

(2) Every triangle attached to a corner is replaced by a refinement, where the 
edges incident to the corner are split as specified by the scaling factor a. This 
process is repeated recursively until the desired nesting level u is reached. 
Note that the mesh may include quadrilaterals after refinement. 

Since the choice of the initial mesh affects strongly the refinement process, it is 
advisable to test with different choices. Naturally, one would want the initial mesh 
to be minimal, that is, having the smallest possible number of elements yet providing 
support for the refinement. This is why initial meshes are sometimes referred to as 
minimal meshes. 

In Figure [3] a challenging example is shown. In this case the large variation of 
the edge lengths is addressed by adding a refinement step to the construction of the 



10 



H.HAKULA, A.RASILA, AND M.VUORINEN 




initial mesh. A detail of the initial mesh is given in Figure |4] along with the final 
mesh. 

3.9. p-Distribution. If so desired, the polynomial degree can be varied over the 
mesh. In Theorem [l] the degrees are connected to the refinement layers. At every 
step in the mesh grading process a new layer of elements is created. Those elements 
belonging to the final layer are assigned the lowest polynomial degree. 

This layer representation can always be recovered by considering graph-distances 
in the graph of the dual mesh, that is, in the graph where the elements are vertices 
and edges indicated a shared edge in the actual mesh. Elements adjacent to points 
of singularity are considered equivalent to those created at the final step. 



4. The case of a rectangle 
The first tests with the hp-FEM software were made for the case of the exterior 



modulus of a rectangle and checked against the Duren-Pfaltzgraff formula (2.5). For 
a convenient parametrization of the computation, the vertices of the rectangle were 
chosen to be the points 1, e**, — 1, — e**, t e (0, 7r/2] of the unit circle. In this case, 
the "interior" modulus of the rectangle is tan(t/2) . It is equal to the modulus of the 
family of curves joining the sides [1, e**] and [—1, — e**] and lying in the interior of the 



rectangle. The formula (2.5) now gives the corresponding exterior modulus as 



2X{k) ' ^ Vtan(t/2) 

For the computation, we carried out the inversion z ^ 1/z in the unit circle 
which keep all the points of the unit circle fixed and transforms the exterior modu- 
lus problem for the rectangle to the "interior" modulus problem of a plane domain 
bounded by four circular arcs, see Figure [5j These circular arcs are the images of 
the sides of the rectangle under the inversion. The results turned out to be quite 
accurate, with a typical relative error of the order 10^^" . 
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Table 1. Exact values of the moduli of e**, — 1, — e**) given by 
(2.5) and the errors of computational results of the hp-method, p = 20, 
the AFEM method and the SC Toolbox. The errors are obtained 



by comparing with the exact formula (2.5). The errors are given as 
Rogio lerror]]. 



k 


exact(t = A;7r/12) 


Error [hpFEM] 


Error [AFEM] 


Error [SCT] 


1 


1.50290233467 


-9 


-6 


-9 


2 


1.31044063554 


-9 


-6 


-9 


3 


1.20035166917 


-9 


-6 


-10 


4 


1.12114255114 


-10 


-6 


-9 


5 


1.05681535228 


-10 


-6 


-13 


6 


1. 


-10 


-6 


-15 



5. Side sliding conjecture 
In this section we will study side sliding problem. 

5.1. The side sliding problem. Consider the problem of finding the minimal 
exterior modulus of the polygonal quadrilateral with vertices 0, 1, a = t + ih,b = 
t — s + ih when h,s > are fixed and t varies. We consider the question of computing 
the modulus of the family F of curves joining the opposite sides [1, a] and [b, 0] outside 
the quadrilateral. Our first step is to reduce the problem to an equivalent problem 
such that three of the points are on the unit circle. Note that this setting is valid 
only if zq is inside the quadrilateral. Indeed, for every choice of h and s this condition 
defines an upper limit for the value of t. 

5.2. Side sliding conjecture. The least value of the exterior modulus is attained 
when t = (1 + s)/2 . For t < (1 + s)/2 the modulus is a decreasing function of t . 

5.3. Numerical experiments on side sliding conjecture. In Figure [7] we show 
a graph of the exterior module as a function of the parameter t G [0.5,2.5], when 
h = l,s = 2. The computation was carried out with SC Toolbox, hp-FEM, and 
AFEM and for the range of computed values, the respective graphs coincide. For the 
SC Toolbox and the hp-FEM the reciprocal estimate for the error was smaller than 
10-*^ and for AFEM 10"^ 



6. The case of a symmetric hexagon 

Suppose that Q{a,b,0, 1) is a quadrilateral in the upper half plane. Then the 
closed polygonal line a, b, 0,6, a, 1, a defines a hexagon H = Q U Q symmetric with 
respect to the real axis. Map the complement of H onto C\{ — 1— t, 1 + t} by a 
conformal map h such that h{0) = —1—t, h{b) = h{b) = —1, h{a) = h(a) = 1, h{l) = 
1 +t where t > depends on the point configuration a, 6, 0, 1 . It is clear by symmetry 
that 



(6.1) 



2M(A+) = M(A) 
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(a) fc = 1 



(b) /c = 2 



{C)k 




(D) k 



(E) k - 



(f) k = 6 



Figure 5. Circular arc domains used for the hp-FEM computations 
of the values in Table [Tj The scale varies from picture to picture. 

where 

A = A([-l-t,-l], [l,l+t];C) and A+ = A([-l-t,-l], [l,l+t];{z : Im^ > 0}) 
Because of the conformal invariance of the modulus we also have 
(6.2) 2M{h~\A+)) = M{h-\A)) . 

Applying this formula with ( 2.5[ ) we see that 

X'{k) ,^ , _i / 1 



(6.3) 



M(r. 



4:X{k) 



k = ip 



2h 
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Figure 6. The circular arc domains of Figure |5| in the same scale. 

1.005 
1 

0.995 

0.99 
0.985 

0.98 
0.975 

0.5 1 1.5 

Figure 7. Side Sliding Conjecture: Dependence of the exterior mod- 
ulus on parameter t with h = l,s = 1 . Maximum is reached at 
t = (1 + s)/2 = 1, as predicted by the conjecture. 

where for h > 0, 

(6.4) r+ = A([0,i/i],[l,l + i/i];C+\[0,l] x [0,/i]). 

This formula can be checked by using the SC Toolbox to construct the above 
conformal mapping h . The tests we carried out for h = 0.2, 0.3, 0, 4 and 0.5 In these 
cases the reciprocal estimate for error was smaller than 10~^ . 
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7. General quadrilateral 

The exterior modulus of the quadrilateral Q with vertices a, 6, c, d is considered 
in this section, i.e., we compute 

I Vul^ (ia; dy 



over the complement of the quadrilateral when u is the solution of the Laplace equa- 
tion in the complement of the quadrilateral with Dirichlet values 1 and on the sides 
[6, c] and [d, a] , respectively, and the Neumann value on the sides [a, h] and [c, d] . 
Here we allow the boundary of the quadrilateral dQ, be a parametrized curve 7(i), 

t e [-1,1]. 

In Figure |8] an overview of the standard FEM approach is given. Using higher- 
order elements one can stretch the domain without introducing significant number 
of elements. Singularities at the corner point must be accounted for in the grading 
of the mesh. Since both the circle and the square cases are symmetric, the exterior 
modulus is exactly 1, and furthermore the potential value at infinity or the far- field 
value is exactly 1/2. 

In Tables [2} [3| and [6] results on two polygonal quadrilaterals 

• Quadrilateral A: {0, 1, (28/25, 69/50), (-19/25, 21/25)}, 

• Quadrilateral B: {0, 1, (42/25, 4), (-3/25, 21/25)}, 

are presented. The exterior modulus has been computed using three methods as an 
equivalent interior modulus problem, and also in truncated domain. In the interior 
case, both SC Toolbox and /ip-FEM give similar results, but AFEM in its standard 
setting does not reach the desired accuracy. This is probably due to the adaptive 
scheme failing in the presence of cusps in the domain. Tables |2] and [3] indicate that 
large exterior angles are the most significant source of errors in the FEM solutions, 
as expected. In the rather benign setting of the Quadrilateral A, SC Toolbox and 
both the internal and external /ip-FEM versions have the same accuracy, but in the 
case of Quadrilateral B, we see gradual loss of accuracy in the FEM solutions. 

Finally, we consider two flower domains, that is, quadrilateral domains with the 
boundary 7(t) = r(t)e** , r{t) = 4/5+(l/5) cos(r;,7rt) and corners at t = —1, —1/4, 0, 1/2. 
For the Quadrilateral C we choose n = A and for D we choose n = 8. These do- 
mains have the useful property that the exterior problem can easily be converted to 
a corresponding interior problem of the domain with boundary l/'y{t). Since these 
domains cannot be handled using the SC Toolbox, we take the interior solution as 
the reference. Tables |4] and |5] show that we can obtain results of high accuracy also 
in traditionally challenging domains. 

It turns out that besides the actual value of the exterior modulus one can also 
determine the value of the far-field potential. Either one can determine the value of 
the potential at the reflection point of the interior problem, i.e., at the origin, or sim- 
ply evaluate the solution of the exterior problem at the farthest point. Remarkably, 
the truncated domain results agree well with the (theoretically) exact results of the 
equivalent inner modulus problems (Table [6]). In Figures |9 12 we show comparisons 



of the interior and exterior potential fields. For the two polygonal quadrilaterals, 
the corresponding contour lines and the location of the origin in the interior case 
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(a) Exterior domain 
with radius > 1000000. 




(b) Zoom of the mesh 
in the case of a circle. 



(c) Zoom of the mesh 
in the case of a square. 




Figure 8. Exterior modulus over the exterior domain. 
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Method 


Exterior Modulus 


Error ( 


2.13 


) 


Relative Error 


SC Toolbox 
AFEM 
hp-FEM (Interior) 
hp-FEM (Exterior) 


0.9923416323 
0.9923500126 
0.9923416332 
0.9923416332 


-9 
-4 
-9 
-9 


-5 
-9 
-9 



Table 2. Quadrilateral A: {0, 1, (28/25, 69/50), (-19/25, 21/25)}. 



The values obtained with SC Toolbox are used as reference. The errors 
are given as [log^^o |error|]. 



Method 


Exterior Modulus 


Error ( 


2.13 


) 


Relative Error 


SC Toolbox 
AFEM 
hp-FEM (Interior) 
hp-FEM (Exterior) 


0.9592571721 
0.9593012739 
0.9592571731 
0.9592572007 


-9 
-4 
-8 
-7 


-4 
-8 
-7 



Table 3. Quadrilateral B: {0, 1, (42/25,4), (-3/25,21/25)}. The val- 
ues obtained with SC Toolbox are used as reference. The errors are 
given as flogj^g |error|]. 



Method 


Exterior Modulus 


Error ( 


2.13 


) 


Relative Error 


hp-FEM (Interior) 
hp-FEM (Exterior) 


0.8196441884805177 
0.8196441926483611 


-14 
-8 


-8 



Table 4. Quadrilateral C: 7(t) = r{t)e'\r{t) = 4/5 + (1/5) cos(47rt) 
and corners at t = —1, —1/4, 0, 1/2. The values obtained with hp-FEM 
(Interior) are used as reference. The errors are given as [log^Q |error|]. 



Method 


Exterior Modulus 


Error ( 


2.13 


) 


Relative Error 


hp-FEM (Interior) 
hp-FEM (Exterior) 


0.9122187602015264 
0.9122187628550672 


-10 
-8 


-8 



Table 5. Quadrilateral D: 7(t) = r{t)e'\r{t) = 4/5 + (1/5) cos(87rt) 
and corners at t = —1, —1/4, 0, 1/2. The values obtained with hp-FEM 
(Interior) are used as reference. The errors are given as [log^^g |error|]. 



are indicated. In the general case, prediction of the far-field value based solely on 
geometric arguments is an open problem. 

We note, that for both Quadrilateral C and D, the interior and exterior capacities 
are equal. This invariance is new and has not been reported in the literature before. 
It is crucial that the four corners are chosen from extremal points, that is, local 
minima and maxima of the radius. 
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Quadrilateral 


hp-FEM (Interior) 


hp-FEM (Exterior) 


Relative Error 


A 


0.5281867366243582 


0.5281867468410989 


-7 


B 


0.6659476737428786 


0.6659476800244547 


-8 


C 


0.5873283399651075 


0.5873283469398137 


-7 


D 


0.5398927341965689 


0.5398927414203410 


-7 



Table 6. Comparison of the computed values of the potential at in- 



finity. The errors are given as [log^g |error|]. 




Figure 9. Quadrilateral A: Correspondence of the potential contours 
between the inner (A) and outer (B) solutions. Shown are the potential 
levels u{z) = 0, 1/10, . . . , 1, and u{0). Corresponding contours can be 
identified by matching the shadings of the regions in between. 



8. Conclusions 

In this study we have shown that three different algorithms, AFEM, SC Tool- 
box and hp-FEM, can all be effectively used for computation of the exterior modulus 
of a bounded quadrilateral. The problem is first reduced to a Dirichlet-Neumann 
problem for the Laplace equation. In our earlier paper [13] we introduced the recip- 
rocal identity as an error estimate for the inner modulus computation and here we 
demonstrate that the same method applies to error estimation for the exterior mod- 
ulus as well. As far as we know, there are very few numerical or theoretical results 
on the exterior modulus in the literature. We compare our numerical results to the 
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(a) Contours of the inner prob- (b) Contours of 

lem. Origin is indicated with a the outer prob- 

dot. lem. Note the 

contours extend- 
ing to infinity. 



Figure 10. Quadrilateral B: Correspondence of the potential con- 
tours between the inner (A) and outer (B) solutions. Shown are the 
potential levels u{z) = 0,1/10,...,!, and u(0). Corresponding con- 
tours can be identified by matching the shadings of the regions in 
between. 




(a) Potential field of (b) Potential field of 

the inner problem. the outer problem. 



Figure 11. Quadrilateral C: The potential field of the inner (A) and 
outer (B) solutions. 

Duren-Pfaltzgraff formula for the exterior modulus of a rectangle and observe that 
our results agree with it, within the limits provided by the reciprocal error estimate. 
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(a) Potential field of 
the inner problem. 



(b) Potential field of 
the outer problem. 



Figure 12. Quadrilateral D: The potential field of the inner (A) and 
outer (B) solutions. 



A problem of independent interest is the value of the potential function at 
infinity. We study this problem for the exterior modulus of a polygonal quadrilateral 
and solve it by mapping the exterior domain onto a bounded domain by inversion 
and then computing the value of the potential function of the corresponding interior 
modulus problem at the image point of the point at infinity. 

We anticipate that the problem of exterior modulus could also be studied using 
integral equation methods. 
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